module mesh_mod

	use const_mod
	use namelist_mod

	implicit none

	private
	public mesh_type
	public mesh

	type mesh_type
		integer lon_halo_width
		integer lat_halo_width
		integer lev_halo_width
		integer num_full_lon
		integer num_half_lon
		integer num_full_lat
		integer num_half_lat
		integer num_full_lev
		integer num_half_lev
		integer full_lon_ibeg
		integer full_lon_iend
		integer full_lat_ibeg
		integer full_lat_iend
		integer full_lev_ibeg
		integer full_lev_iend
		integer half_lon_ibeg
		integer half_lon_iend
		integer half_lat_ibeg
		integer half_lat_iend
		integer half_lev_ibeg
		integer half_lev_iend
		integer full_lon_lb
		integer full_lon_ub
		integer half_lon_lb
		integer half_lon_ub
		integer full_lat_lb
		integer full_lat_ub
		integer half_lat_lb
		integer half_lat_ub
		integer full_lev_lb
		integer full_lev_ub
		integer half_lev_lb
		integer half_lev_ub
		real(r8) start_lon
		real(r8) end_lon
		real(r8) start_lat
		real(r8) end_lat
		real(r8) dlon
		real(r8) dlat
		real(r8), allocatable, dimension(:) :: full_lon
		real(r8), allocatable, dimension(:) :: half_lon
		real(r8), allocatable, dimension(:) :: full_lat
		real(r8), allocatable, dimension(:) :: half_lat
		real(r8), allocatable, dimension(:) :: full_lev
		real(r8), allocatable, dimension(:) :: half_lev
		real(r8), allocatable, dimension(:) :: full_dlev
		real(r8), allocatable, dimension(:) :: half_dlev
		real(r8), allocatable, dimension(:) :: full_lon_deg
		real(r8), allocatable, dimension(:) :: half_lon_deg
		real(r8), allocatable, dimension(:) :: full_lat_deg
		real(r8), allocatable, dimension(:) :: half_lat_deg
		real(r8), allocatable, dimension(:) :: full_cos_lon
		real(r8), allocatable, dimension(:) :: half_cos_lon
		real(r8), allocatable, dimension(:) :: full_sin_lon
		real(r8), allocatable, dimension(:) :: half_sin_lon
		real(r8), allocatable, dimension(:) :: full_cos_lat
		real(r8), allocatable, dimension(:) :: half_cos_lat
		real(r8), allocatable, dimension(:) :: full_sin_lat
		real(r8), allocatable, dimension(:) :: half_sin_lat
		! Edge length
		real(r8), allocatable, dimension(:) :: de_lon
		real(r8), allocatable, dimension(:) :: de_lat
		real(r8), allocatable, dimension(:) :: le_lat
		real(r8), allocatable, dimension(:) :: le_lon
		! Cell area
		real(r8), allocatable, dimension(:) :: area_cell
	contains
	  procedure :: init  => mesh_init
	  procedure :: clear => mesh_clear
	  final :: mesh_final
	end type mesh_type

	type(mesh_type), target :: mesh

contains
	
	subroutine mesh_init(this, lon_halo_width, lat_halo_width, lev_halo_width)

		class(mesh_type), intent(inout) :: this
		integer, intent(in), optional :: lon_halo_width, lat_halo_width, lev_halo_width
		integer i, j

		call this%clear()

		this%num_full_lon = num_lon
		this%num_half_lon = num_lon 
		this%full_lon_ibeg = 1
		this%full_lon_iend = this%num_full_lon
		this%half_lon_ibeg = 1
		this%half_lon_iend = this%num_half_lon
		this%num_full_lat  = num_lat
		this%num_half_lat  = num_lat - 1
		this%full_lat_ibeg = 1
		this%full_lat_iend = this%num_full_lat
		this%half_lat_ibeg = 1
		this%half_lat_iend = this%num_half_lat
		this%num_full_lev  = num_lev
		this%num_half_lev  = this%num_full_lev + 1
		this%full_lev_ibeg = 1
		this%full_lev_iend = this%num_full_lev
		this%half_lev_ibeg = 1
		this%half_lev_iend = this%num_half_lev

		this%lon_halo_width = merge(lon_halo_width, 1, present(lon_halo_width))
		this%lat_halo_width = merge(lat_halo_width, 1, present(lat_halo_width))
		this%lev_halo_width = merge(lev_halo_width, 2, present(lev_halo_width))

		this%start_lon      = 0.0_r8
		this%end_lon        = pi2
		this%start_lat      = -pi05
		this%end_lat        = pi05

	  this%dlon           = (this%end_lon - this%start_lon) / this%num_half_lon
	  this%dlat           = (this%end_lat - this%start_lat) / this%num_half_lat

	  this%full_lon_lb    = this%full_lon_ibeg - this%lon_halo_width
	  this%full_lon_ub    = this%full_lon_iend + this%lon_halo_width
	  this%full_lat_lb    = this%full_lat_ibeg - this%lat_halo_width
	  this%full_lat_ub    = this%full_lat_iend + this%lat_halo_width
	  this%half_lon_lb    = this%half_lon_ibeg - this%lon_halo_width
	  this%half_lon_ub    = this%half_lon_iend + this%lon_halo_width
	  this%half_lat_lb    = this%half_lat_ibeg - this%lat_halo_width
	  this%half_lat_ub    = this%half_lat_iend + this%lat_halo_width
		this%full_lev_lb    = this%full_lev_ibeg - this%lev_halo_width
		this%full_lev_ub    = this%full_lev_iend + this%lev_halo_width
		this%half_lev_lb    = this%half_lev_ibeg - this%lev_halo_width
		this%half_lev_ub    = this%half_lev_iend + this%lev_halo_width

	  allocate(this%full_lon    (this%full_lon_lb:this%full_lon_ub))
	  allocate(this%half_lon    (this%half_lon_lb:this%half_lon_ub))
	  allocate(this%full_lat    (this%full_lat_lb:this%full_lat_ub))
	  allocate(this%half_lat    (this%half_lat_lb:this%half_lat_ub))
	  allocate(this%full_lev    (this%full_lev_ibeg:this%full_lev_iend))
	  allocate(this%full_dlev   (this%full_lev_ibeg:this%full_lev_iend))
	  allocate(this%half_lev    (this%half_lev_ibeg:this%half_lev_iend))
	  allocate(this%half_dlev   (this%half_lev_ibeg:this%half_lev_iend))
	  allocate(this%full_lon_deg(this%full_lon_lb:this%full_lon_ub))
	  allocate(this%half_lon_deg(this%half_lon_lb:this%half_lon_ub))
	  allocate(this%full_lat_deg(this%full_lat_lb:this%full_lat_ub))
	  allocate(this%half_lat_deg(this%half_lat_lb:this%half_lat_ub))
	  allocate(this%full_cos_lon(this%full_lon_lb:this%full_lon_ub))
	  allocate(this%half_cos_lon(this%half_lon_lb:this%half_lon_ub))
	  allocate(this%full_sin_lon(this%full_lon_lb:this%full_lon_ub))
	  allocate(this%half_sin_lon(this%half_lon_lb:this%half_lon_ub))
	  allocate(this%full_cos_lat(this%full_lat_lb:this%full_lat_ub))
	  allocate(this%half_cos_lat(this%half_lat_lb:this%half_lat_ub))
	  allocate(this%full_sin_lat(this%full_lat_lb:this%full_lat_ub))
	  allocate(this%half_sin_lat(this%half_lat_lb:this%half_lat_ub))
	  allocate(this%de_lon      (this%full_lat_lb:this%full_lat_ub))
	  allocate(this%le_lon      (this%full_lat_lb:this%full_lat_ub))
	  allocate(this%de_lat      (this%half_lat_lb:this%half_lat_ub))
	  allocate(this%le_lat      (this%half_lat_lb:this%half_lat_ub))
	  allocate(this%area_cell   (this%full_lat_lb:this%full_lat_ub))

	  do i = this%full_lon_lb, this%full_lon_ub
	  	this%full_lon(i) = this%start_lon + (i - 1) * this%dlon
	  end do
	  
	  do i = this%half_lon_lb, this%half_lon_ub
	  	this%half_lon(i) = this%full_lon(i) + 0.5_r8 * this%dlon
	  end do

	  do j = this%half_lat_lb, this%half_lat_ub
	  	this%half_lat(j) = this%start_lat + (j - 0.5) * this%dlat
	  end do

	  do j = this%full_lat_lb, this%full_lat_ub - 1
	  	this%full_lat(j) = this%half_lat(j) - 0.5_r8 * this%dlat
	  end do
	  this%full_lat(this%full_lat_ub) = this%half_lat(this%half_lat_ub) + 0.5_r8 * this%dlat

	  ! this%full_lat(1) = this%start_lat
	  ! do j = 2, this%num_full_lat - 1
	  ! 	this%full_lat(j) = this%full_lat(j-1) + this%dlat
	  ! end do
	  ! this%full_lat(this%num_full_lat) = this%end_lat

	  do i = this%full_lon_lb, this%full_lon_ub
	  	this%full_sin_lon(i) = sin(this%full_lon(i))
	  	this%full_cos_lon(i) = sin(this%full_lon(i))
	  	this%full_lon_deg(i) = this%full_lon(i) * deg
	  end do

	  do i = this%half_lon_lb, this%half_lon_ub
	  	this%half_cos_lon(i) = cos(this%half_lon(i))
	  	this%half_sin_lon(i) = sin(this%half_lon(i))
	  	this%half_lon_deg(i) = this%half_lon(i) * deg
	  end do

	  do j = this%full_lat_lb, this%full_lat_ub
	  	! if (this%full_lat(j) >= -pi05 .and. this%full_lat(j) <= pi05) then
	  		this%full_cos_lat(j) = cos(this%full_lat(j))
	  		this%full_sin_lat(j) = sin(this%full_lat(j))
	  		this%full_lat_deg(j) = this%full_lat(j) * deg
	  	! end if
	  end do

	  do j = this%half_lat_lb, this%half_lat_ub
	  	if (this%half_lat(j) >= -pi05 .and. this%half_lat(j) <= pi05) then
	  		this%half_cos_lat(j) = cos(this%half_lat(j))
	  		this%half_sin_lat(j) = sin(this%half_lat(j))
	  		this%half_lat_deg(j) = this%half_lat(j) * deg
	  	end if
	  end do

	  do j = this%full_lat_ibeg + 1, this%full_lat_iend - 1
	  	this%de_lon(j) = radius * this%full_cos_lat(j) * this%dlon
	  	this%le_lon(j) = radius * this%dlat
	  end do

	  this%de_lon(this%full_lat_ibeg) = 0.0_r8
	  this%le_lon(this%full_lat_ibeg) = 0.0_r8
	  this%de_lon(this%full_lat_iend) = 0.0_r8
	  this%le_lon(this%full_lat_iend) = 0.0_r8

	  do j = this%half_lat_ibeg, this%half_lat_iend
	  	this%de_lat(j) = radius * this%dlat
	  	this%le_lat(j) = radius * this%half_cos_lat(j) * this%dlon
	  end do

	  do j = this%full_lat_ibeg + 1, this%full_lat_iend - 1
	  	this%area_cell(j) = radius**2 * this%full_cos_lat(j) * this%dlon * this%dlat
	  end do

	  this%area_cell(this%full_lat_ibeg) = radius**2 * this%dlon * (this%half_sin_lat(this%half_lat_ibeg) + 1.0_r8)
	  this%area_cell(this%full_lat_iend) = radius**2 * this%dlon * (1.0_r8 - this%half_sin_lat(this%half_lat_iend))
	
	end subroutine mesh_init

	subroutine mesh_clear(this)

		class(mesh_type), intent(inout) :: this

		if (allocated(this%full_lon        )) deallocate(this%full_lon    )
		if (allocated(this%half_lon        )) deallocate(this%half_lon    )
		if (allocated(this%full_lat        )) deallocate(this%full_lat    )
		if (allocated(this%half_lat        )) deallocate(this%half_lat    )
		if (allocated(this%full_lev        )) deallocate(this%full_lev    )
		if (allocated(this%full_dlev       )) deallocate(this%full_dlev   )
		if (allocated(this%half_lev        )) deallocate(this%half_lev    )
		if (allocated(this%half_dlev       )) deallocate(this%half_dlev   )
		if (allocated(this%full_lon_deg    )) deallocate(this%full_lon_deg)
		if (allocated(this%half_lon_deg    )) deallocate(this%half_lon_deg)
		if (allocated(this%full_lat_deg    )) deallocate(this%full_lat_deg)
		if (allocated(this%half_lat_deg    )) deallocate(this%half_lat_deg)
		if (allocated(this%full_cos_lon    )) deallocate(this%full_cos_lon)
		if (allocated(this%full_sin_lon    )) deallocate(this%full_sin_lon)
		if (allocated(this%half_cos_lon    )) deallocate(this%half_cos_lon)
		if (allocated(this%half_sin_lon    )) deallocate(this%half_sin_lon)
		if (allocated(this%full_cos_lat    )) deallocate(this%full_cos_lat)
		if (allocated(this%full_sin_lat    )) deallocate(this%full_sin_lat)
		if (allocated(this%half_cos_lat    )) deallocate(this%half_cos_lat)
		if (allocated(this%half_sin_lat    )) deallocate(this%half_sin_lat)
		if (allocated(this%de_lon          )) deallocate(this%de_lon      )
		if (allocated(this%le_lon          )) deallocate(this%le_lon      )
		if (allocated(this%de_lat          )) deallocate(this%de_lat      )
		if (allocated(this%le_lat          )) deallocate(this%le_lat      )
		if (allocated(this%area_cell       )) deallocate(this%area_cell   )

	end subroutine mesh_clear

	subroutine mesh_final(this)

		type(mesh_type), intent(inout) :: this

		call this%clear()

	end subroutine mesh_final

end module mesh_mod